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' Abstract 

(N 

^ A discontinuous Galerkin method for approximating the Vlasov-Poisson system of 

equations describing the time evolution of a coUisionless plasma is proposed. The 
method is mass conservative and, in the case that piccewise constant functions are 
CO used as a basis, the method preserves the positivity of the electron distribution 

function and weakly enforces continuity of the electric field through mesh interfaces 
and boundary conditions. The performance of the method is investigated by com- 
puting several examples and error estimates associated system's approximation are 

•• stated. In particular, computed results are benchmarked against established theoret- 

• ical results for linear advection and the phenomenon of linear Landau damping for 

both the Maxwell and Lorcntz distributions. Moreover, two nonlinear problems are 

^ considered: nonlinear Landau damping and a version of the two-stream instability 

are computed. For the latter, fine scale details of the resulting long-time BGK-like 
state are presented. Conservation laws are examined and various comparisons to 
theory are made. The results obtained demonstrate that the discontinuous Galerkin 
method is a viable option for integrating the Vlasov-Poisson system. 
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1 Introduction 



The single species Vlasov-Poisson system is a nonlinear kinetic system that 
models time evolution of a collisionless plasma consisting of electrons and a 
uniform background of fixed ions under the effects of a self-consistent elec- 
trostatic field and possibly an externally supplied field. The Vlasov equation 
models the transport of the electrons that are coupled to the electrostatic 
potential through Poisson's equation. 

The unknown electron distribution function, a phase space density, is denoted 
by / = fix, V, t), where the independent variables x, v, t are position, velocity, 
and time, respectively. For a given t, the quantity f{x,v,t)dxdv denotes the 
number of electrons contained in the infinitesimal phase space volume element 
dxdv centered about {x,v) at time t. Upon a proper renormalization, / can 
be interpreted as a probability distribution function for the electrons over the 
phase space. 

The Vlasov-Poisson system has solutions that can exhibit a variety of dynam- 
ical phenomena [691I7E1] . One of the most well-known effects is filamentation 
or 'phase mixing' as it is sometimes called, which occurs when different char- 
acteristics surfaces associated to the nonlinear transport (Vlasov) equation 
wrap in the phase space. This effect results in stiff gradients, since / generally 
takes on disparate values along different characteristics. 

Related to filamentation is the interesting phenomenon of Landau damping 
^2] : electrostatic disturbances can be interpreted as the interaction between 
plasma waves and electrons with a resulting net energy transfer from the wave 
to electrons leading to an exponential collisionless damping of the electric field 
modes in time. Because of such phenomena, the Vlasov-Poisson equation can 
be challenging to simulate numerically. 

Existing numerical techniques for solving the Vlasov-Poisson system can be 
divided into two groups: (i) those that approximate the system in the phase 
space directly and (ii) those that transform the system into a different coordi- 
nate space. The numerical approaches that treat the phase space directly do 
not, however, usually involve discretizing the Vlasov equation. Rather most 
of these methods take advantage of the characteristic structure of the Vlasov 
equation, which implies that the plasma particles evolve along trajectories that 
satisfy a given set of ordinary differential equations. The most widely used par- 
ticle scheme is the Particle- in-Cell (PIC) method [9)1381)2 7j . which represents 
ensembles micro-particles as a finite number of macro-particles. Each macro- 
particle is then assumed to evolve along a characteristic trajectory, where the 
electric field defining the trajectory is computed via any standard scheme. 
The PIC method seems to give reasonable results in cases where the tail of 
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the distribution is negligible and a large number of particles is not necessary. 
Otherwise, the method suffers from numerical noise that is proportional to 
I/a/ti, where n is the number of particles. 

Other methods based on the discretization of the phase space have also been 
proposed. In [11], an operator splitting method was introduced and shown 
to be both efficient and accurate for solving a wide range of problems. A 
continuous finite element method was developed in [ 72|75] and was shown to 
achieve results similar to those obtained in pJ^J . A positive and mass conserva- 
tive scheme was employed in |33] to solve both linear and nonlinear damping 
problems in one- and two-dimensional physical space. This method is defined 
at a given time step by first building a piecewise constant approximation over 
a mesh of the phase space using the approximation obtained from the previ- 
ous time step and two correction terms whose values are found by solving two 
fixed point problems. The piecewise approximation is then used in conjunction 
with a slope limiter to reconstruct a local polynomial approximation of / for 
each cell in the mesh. 

Transform methods based on Fourier or Hermite series have also been used, 
e.g. in [36f3ll29] and more recently in |3T]. In |^OKT] . the phase space was 
transformed using a Fourier transform and a splitting method was employed to 
advance the approximation in time. This method also included a filamentation 
filtering step for the purpose of smoothening the filaments. The numerical 
results obtained using this method seem reasonable only for problems in which 
filamentation is not a dominant effect, where perhaps the nonlinearity either 
slows down or prevents the onset of filamentation. However, this method may 
be inadequate for problems where the physics of interest depends upon the 
filamentary nature of the distribution, such as is the case for Landau damping. 

The objective of this paper is to propose a coupled Upwind Penalty Galerkin 
(UPG) method for the approximation of the Vlasov-Poisson system and to 
evaluate its numerical efficacy. Our UPG method gives a unified approach for 
approximating both the hyperbolic and elliptic parts of the Vlasov-Poisson 
system. Specifically, the Vlasov equation is discretized using the standard Up- 
wind Galerkin (UG) scheme for conservation laws [Mll23ll22j and the Poisson 
equation is discretized using one of the three Discontinuous Galerkin (DG) 
interior penalty available schemes [70|59II63] . Stability and convergence esti- 
mates for the UPG method were presented in [37] for the six-dimensional 
phase space In the same reference, the method was shown to be both locally 
and globally mass conservative. 

More specifically, the semi-discrete UPG approximation fh to the electron 
distribution function is defined to be the solution of a first-order, nonlinear, 
ordinary differential equation (ODE) system. Moreover, it has been shown that 
the method preserves positivity of fh when piecewise constants basis functions 
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are used to approximate the solution to the Vlasov-Poisson system. 

In this manuscript we show the numerical efficacy of the DG method by 
performing accuracy, convergence, and conservation tests on computed UPG 
approximate solutions to a variety of linear and nonlinear problems where 
sufficiently data or information of the true nonlinear solutions has been es- 
tablished. The computed results for these problems are benchmarked against 
known theoretical results and are compared to results obtained using estab- 
lished methods. 

For computing plasma problems the UPG method offers several advantages. 
In particular, the local nature of the method facilitates adaptive mesh refine- 
ments with easy adaptation to parallelization techniques. By taking advantage 
of these beneffis, regions of the phase space where the electron distribution 
experiences strong filamentation or boundary layer effects can be resolved by 
local mesh refinements. The discontinuous nature of the method also helps 
to resolve the stiff gradients associated with filamentation, since requiring the 
approximation to be continuous in these cases can be too restrictive and typ- 
ically lead to excessive numerical diffusion and oscillatory behavior. Due to 
the fact that the method imposes boundary conditions weakly, a variety of 
boundary conditions can easily be accommodated. 

Recently, an alternative DG formulation for the Boltzmann-Poisson system 
was introduced in [T^13lll4jT^6j for simulations of hot electron transport for 
one and two space dimensional nanometer scale devices, where kinetic correc- 
tions are known to be very significant. In [20] a DG scheme was constructed for 
the Vlasov-Boltzmann equation by means of a maximum principle satisfying a 
limiter for conservation laws [TlfTS] , and it was shown to be high order accu- 
rate and positivity-preserving, not only for piecewise constant basis functions 
but also for higher order polynomial approximations. Consequently a new ap- 
proximation of the Vlasov system based on a UPG scheme for higher order 
polynomial basis functions is currently being developed and tested (HUT]. 
Thus, the DG method is well-suited to approximate a range of plasmas span- 
ning from the coUisionless to the highly collisional regimes. 

For completeness of this introduction we include some historical remarks re- 
garding DG methods. The initial development of these numerical approximat- 
ing schemes for hyperbolic and elliptic equations occurred independently, but 
nearly simultaneously. In 1973, the ffist DG scheme for linear hyperbolic equa- 
tions was introduced by Reed and Hill for approximating a neutron transport 
equation [57] . This work was followed by Lesaint and Raviart j33] in 1974, 
where a priori error estimates were proved for the DG method applied to 
two-dimensional, linear hyperbolic problems. The DG schemes for hyperbolic 
problems were further studied in the series of papers [Mll23f22] . which culmi- 
nated in the introduction of the local discontinuous Galerkin (LDG) method 
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|25j . The generality of the LDG method was further extended to the multi- 
dimensional setting under more relaxed assumptions on the data in ^T\. One 
of the first DG schemes for approximating the solutions to second-order ellip- 
tic equations was introduced in 1971 by Nitsche [55] where Dirichlet bound- 
ary conditions were enforced weakly rather than strongly through the use 
of a penalty term. Shortly thereafter, applications of the penalty method to 
Laplace's equation were proposed by Babuska et al. in [H[5|6] . 

The use of penalty terms across interior faces as a means of enforcing continu- 
ity among adjacent elements was introduced in [7D] and [7T] using a symmetric 
interior penalty Galerkin (SIPG) finite element method. A non-symmetric in- 
terior penalty method (NIPG) similar to the SIPG method was proposed and 
analyzed in |59j. The incomplete interior penalty method (IIPG) was intro- 
duced in j63ii28p64j and is very similar to the SIPG and NIPG methods. 

The outline of this paper is as follows. In Section [2| we describe the Vlasov- 
Poisson system and give a brief discussion of Landau damping. In Section 
[3} the UPG method for the approximation of the Vlasov-Poisson system is 
derived and the error estimates associated to the approximation of the system 
are stated. In Section |4| several numerical simulation results are presented and 
analyzed, including a study of the free streaming operator (simple advection). 
Landau damping for Maxwellian and Lorenztian equilibria, strong nonlinear 
Landau damping and a careful study of a symmetric two-stream instability, 
all for the case of repulsive potential forces. In Section |5j we comment on the 
efficiency of the UPG method, draw some conclusions, and remark on future 
work. Finally, an Appendix dedicated to the analysis of dispersion relations 
for the Lorentzian and two-stream equilibria is included. 



2 The Vlasov-Poisson System 

The Vlasov-Poisson system considered in this work has been scaled by the 
usual characteristic time and length scales, i.e., time is scaled by the inverse 
plasma frequency u~^, length by the Debye length Ad, and velocity accordingly 
by a thermal velocity Vth = ujp^D- 

Using this nondimensionalization, the Vlasov-Poisson problem is as follows: 
For the divergence free vector in {x, v) phase space 

a{x,v,t) = {v,-E{v,t)) for (x, t;, t) G ^] X [0, T] , (1) 

find the electron distribution function f{x,v,t) and the electric field E{x,t) 
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with corresponding electrostatic potential t) such that, for fixed T > 0, 

= dtf + div{af) = dJ + vVJ-E-VJ-, infix (0,T], (2a) 
£;=-V,$; -A,$ = l-/ fdv, infi,x(0,r], (2b) 

subject to an initial condition on / and boundary conditions on / and The 
domain of definition of the initial boundary value problem Q :— fi-^ x R", 
where the physical domain C M"^ can either be bounded or all of M"^ with 
n = 1,2,3. The boundary condition given for / depends on Q^- If = 
then the condition / ^ must be enforced both as \x\ oo and as \v\ oo. 
If Jlj. is bounded, then a condition must be imposed on / along the infiow 
boundary F^, defined by 

r, = {{x, v) edn:,xW\viy^<0}, (3) 

with being the unit outward normal vector to dflx- Often is given in 
some parts of the boundary and may be periodic in other parts of the boundary 
region T^. In this manuscript wc assume periodic boundary conditions in space 
and the decaying boundary condition in velocity. 

The Poisson equation must also be endowed with spatial boundary conditions, 
either Dirichlet, Neumann, Robin, or periodic, on different disjoint regions of 
the boundary dflx- We denote the Dirichlet portion of the boundary by dflx,D- 
If the measure of the Dirichlet boundary is zero, i.e. |5fia;,D| = 0, and one has 
homogeneous or periodic boundary conditions such that Jqq^ V$ • — 0, 
then in order to maintain a well-posed problem that keeps the existence and 
uniqueness of the corresponding Poisson boundary value problem, one needs 
to add to the solution space the compatibility (neutrality) condition /^j (1 — 
/r" fdv)dx = 0, or equivalently /^^ /r« fdxdv — \^x\ on each connected 
component of the spatial domain fi. 

Macroscopic fluid quantities of interest are easily computed from /. The elec- 
tron density p = p{x,t), current density j = j{x,t), kinetic energy density 
£k{x,t) and electrostatic energy density Se{x,t) are defined by 



p{x,t) 


= j f{x,v,t)dv, 

R" 


(4) 


j{x,t) 


= J vf{x,v,t)dv, 

R" 


(5) 


£k{x,t) 


^ 2 / 1^1^(2^,^^,^)^?^^, 

R" 


(6) 


£e{x,t) 


- l\E{x,t)\'. 


(7) 



These quantities satisfy a number of respective conservation laws (see e.g. 
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|58j). In particular, it is well-known that the Vlasov-Poisson system conserves 
total particle number, momentum, energy, and the Casimir invariants, which 
are given, respectively, by 



N 
P 
H 

C 



f{x,v,t)dxdv= / p{x,t)dx, 
vf{x,v,t)dxdv= / j{x,t)dx, 

If 1 

- / \v\'^ f{x,v,t) dx dv + - 
2 Jq^xR" 2 

(£k{x,t) +£e{x,t)'j dx , 



E{x,t)\^dx 



dx dv 



(9) 
(10) 

(11) 



The notation C{f) in ( [lT| ) refers to an arbitrary function of / and includes the 
'enstrophy' when C{f) = entropy when C{f) = —/In/, or particle number, 
as in (Is]), when C(/) = /. We will check the invariance of these quantities in 



our nonlinear computations, particularly in Section 4.2.2 



Interesting properties of the Vlasov-Poisson system result by considering a 
linear perturbation 6f{x,v,t) to an equilibrium distribution feq{v) over the 
2D-phase space [0,L] x (— oo,cxd), L > 0. Specifically, suppose that / = 
feq + ^f, where feq is a given equilibrium probability distribution, 6f and $ 
are L-periodic in x, and the initial average value of 6f over Q is zero. Equations 



(2a)- (2b) imply that 5f satisfies 



dt{6f)+v{6f),-E{6f), 

E 



eq 



[0,L] X (-00,00) X (0,T] 
[0,L] X (0,T] 

6fdv [0,L]x(0,T] 



Supposing \E{Sf)v\ « 1 and dropping this term from (12a) leads to 



dt{Sf)+v{6f), 
E 



EfL 



eq 



5f dv 



[0,L] X (-00,00) X (o,r] , 
[0,L] X (0,T], 

[0,L] X (0,T]. 



(12a) 
(12b) 

(12c) 



(13a) 
(13b) 

(13c) 



The linear system of (13a)-(13c) was analyzed by using the Laplace transform 



in the famous paper of Landau |42j , by expansion in terms of continuum eigen- 
functions in [69], and by a tailored integral transform introduced in [5T|52] . 
Landau showed that an electric field mode Ek{u,t) decays exponentially in 
the long-time limit, which we investigate in Section [4.1.2 for two well-known 
equilibria: the Maxwellian 



/m = (2vrr) 



^n/2 g-»)2/2T ^ 



(14) 
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and the Lorentzian, 



fLiv)= (15) 



3 Method Formulation 



In this section, we derive the UPG method for the Vlasov-Poisson system. The 
derivation proceeds by first discretizing the Vlasov equation using the stan- 
dard UG discretization for transport equations [21]. Here it is assumed that 
the electric field is given and hence the divergence-free flow field a{x,v,t) = 
E(x, t)), defined in ([T]) for the Vlasov equation (2a), is known. Afterwards, 
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the DG discretization for the Poisson equation is considered. 

It is assumed that any mesh for the phase space, where by mesh we mean a par- 
titioning of the phase space into convex sets called elements, is the Cartesian 
cross product of a mesh for the physical domain and a mesh for the velocity 
domain. Under this assumption, the physical and velocity domains can be in- 
dependently refined. Given a mesh for the phase space, the UG method then 
defines an approximate solution to the true Vlasov solution in such a way that 
at any given time the approximate solution restricted to each element of the 
mesh is a polynomial function. However, the approximate solution is not re- 
quired to be continuous across the intersections of any two adjacent elements, 
so it is a piecewise defined polynomial function with respect to the mesh at 
any given time. 

In order to compute the approximation to the electrostatic potential from 



Poisson's equation (2b), for a given distribution function, one may use one 
of three interior penalty methods that weakly enforce both approximate con- 
tinuity across the interior mesh faces and Dirichlet boundary regions. These 
three alternative methods, symmetric interior penalty Galerkin (SIPG) [TCTfTT] . 
non-symmetric interior penalty Galerkin (NIPG) [SS], and incomplete interior 
penalty Galerkin (IIPG) [63m64] are discussed in detail and a priori error 
estimates for each of them are given in the respective references. The only dif- 
ference among the three methods is in the value of one specific parameter that 
arises in the weak formulation that is common to each of them. Thus, for a 
given space charge function, each penalty Galerkin method defines a piecewise 
polynomial approximation of the true solution to the Poisson equation using 
a mesh in the spatial domain. 

The spatial domain mesh used in the discretization of Poisson's equation is 
required to be the same as that used in the UG discretization of the Vlasov 
equation. This requirement is a practical one, both in terms of analysis and 
implementation. However, the polynomial degree of the potential approxima- 
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tion on a given element of the spatial mesh is not required to equal the degree, 
with respect to x, of the polynomial approximation of the distribution /. 

Thus, the UPG method of approximation to the Vlasov-Poisson system is de- 
fined by coupling together the UG method of approximation to the Vlasov 
equation with interior penalty methods of approximation to the Poisson equa- 
tion. This nonlinear semi-discrete approximation results in a first-order non- 
linear ODE system, the solution of which determines the approximation fh of 
/. The resulting ODE system is readily solved using an explicit conservative 
time-integrator such as the Runge-Kutta method. Moreover, in the process of 
computing fh, both the approximation Eh of the electric field E and the ap- 
proximation $/j of the potential $ are computed by one of the penalty methods 
for the elliptic equation. 



3. 1 Preliminaries 



We assume the computational domain in velocity space is a bounded set VL^ 
and that the approximate solution to f{x,v,t) is assumed to vanish in dQ^. 
It is then, implicitly for our simulation, assumed that the velocity support 
of the approximation to the true solution / of the Vlasov-Poisson system is 
contained in for all times. This is a reasonable assumption for problems with 
spatial periodic boundary conditions, as it is expected that most of the density 
associated with the approximation to / will be contained in a sufficiently large 
fixed set Qy. The error due to this assumption only depends on the density of 
the true solution / computed on the complement of fl^ in M". 



Further, the conservative nature of the transport equation (2a) with the space 
dependent divergence-free fiow field a, motivates us to choose the computa- 
tional scheme as follows: 

Let {Th^}hj:>o be a sequence of successively refined meshes of the bounded 
domain fl^ C M" and let {Th^}h^>o be a sequence of successively refined meshes 
of the bounded domain Q^j C MJ^, where n = 1, 2, 3. Given the meshes Th^ = 
{i^j^l^t^i and 7/i„ = {i^'j^l^t^i, the elements Kj^ and Kj^ comprising each 
of the respective meshes are sets of the following types: intervals, if = 1; 
triangles or quadrilaterals, if = 2; and tetrahedra, prisms or hexahedra, if 
n = 3. The corresponding spatial refinement level and velocity refinement 
level hv are defined by hx = maxj^ {diam{Kj^)} and = max^^ {diam(iC,-^)}, 
respectively. 

A sequence of successively refined meshes {Th}h>o of the, now, computational 
domain f2 = i^a, x f2„ is generated by defining each mesh % = {Kj}^J^i to 
he Th = 7h^ X Th^, where the refinement level is h = [h^^h^). Thus, for any 
given element Kj G Th there exists a unique pair of elements Kj^ G Th^ and 
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Kj^ G Th^ such that Kj = Kj^ x Kj^, which is equivalent to the existence 
of an invertible mapping from j G {1, . . . , Nh} to {jx,jv) ^ {!)••• ; ^ 
{1, . . . , NhJ, where A^^, = N^^^Nh.^. 

The derivation of the following UPG method requires the use of the broken 
Sobolev space H'^lTh}, s > 1/2, which is defined as follows: 

H'{%} = {we L\n) I w^^^ G H'{K,], \JK, e %} , (16) 

i.e., H'^{Th} is the space of those functions that have elementwise weak deriva- 
tives up to, and including, the order s. Then for nonnegative integers and 
r^, the discontinuous approximation space D^'^^^iTh) C H'^{Th} is given by 

D^-'^iTH) = {we H-'{rh} I G Q'^{K,J X Q^^{K,J, ^ K, G } , 

(17) 

where Q^(ii") denotes the space of polynomials on a set K with degree less 
than or equal to r in each variable. Thus, F^{K) C Q''(K), where F'^{K) 
denotes the space of polynomials satisfying that the sum of the degrees of all 
the variables is less than or equal to r. 

The choice of Q^{K) for basis functions is suitable for Cartesian meshes in 
both x-space and w-space, respectively, where trace and inverse inequalities 
that are derived by mapping to the reference element in the approximating 
framework are possible, as first introduced in [37] . 

However, one may use triangles in two-dimensions, and prisms or hexahedra in 
three-dimensions, for both x-space and f-space, for which the natural choice 
of polynomial space would be F'^^{Kj^} x F^''{Kj^}. We point out that these 
selections of approximating spaces is consistent with the divergence free, linear, 
and conservation form of the Vlasov equation, and the fact that the choice 
of the mesh associated with the computational domain Q = fl^ x is a set 
of product mesh elements Kj = Kj^ x Kj,^, for j G {1, . . . , Nh}, {jx,jv) ^ 
{1, . . . , NhJ X {1, . . . , NhJ and A^,, = N^^Nh,. Such is clearly preserved by 
the Vlasov flow and the corresponding approximation results for F^{K) will 
be valid. This choice may be preferable in higher dimensions since the number 
of degrees of freedom of the basis functions in F'^{K) is rn + 1, while for <[Y{K) 
is (r + 1)", for n-dimensional calculations. 

The discontinuous nature of the space H'^{Th} needs the introduction of 
mesh faces. If Kj is a boundary element, then = dKj fl dfl is called 
a boundary mesh face. If Ki and K2 are two intersecting elements whose 
common intersection lies in the interior of Q, then Fk = dKi fl 8X2 is said 
to be an interior mesh face. The set of all mesh faces is denoted by J^h = 
{Fl, . . . , Fp^, -Fp^+i, . . . , Fm^}, where F^ is an interior face ii I < k < Ph and 
a boundary face if -P/i + 1 < A; < M^. Each face Fk G J-'h is associated with a 
unit normal vector z/^. For k > Ph, Vk is chosen to be the outward unit normal 
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to dn. For 1 < A; < P/i, we fix z/fc to be one of tlie two unit normal vectors to 
Fk. For every interior face F^, the elements Ki and K2 will always be used to 
denote the two unique elements such that F^ = KinK2. Moreover, it is always 
assumed that Ki satisfies z/^^ = Uk on F^, where i/^^ denotes the outward unit 
normal to dKi. Then = —Uk on Fk- 

The fact that each mesh 7h is the product of a spatial mesh Th^ and a velocity 
mesh 7h„ gives a specific structure to the boundaries of the elements and to the 
set of mesh faces Th- It follows that for Kj = Kj^ x Kj^ we have dKj = {dKj^ x 
Kj^) U {Kj^ X dKjJ). We denote the set of mesh faces for Th^ and Th^ by J^h^ = 
{Ff , . . . , F^^^ , i^p^^+i, . . . , i^Mft, } and J^h^ = {F^, . . . , F^^^ , i^p^^+i, • • • , F^^J, 
respectively, where F^^ is an interior face if 1 < fc^; < Ph^ and a boundary face 
if Ph^ + 1 < < ^/lo:! ciiid is an interior face if 1 < fc^, < Ph^ and a 
boundary face if Ph^ + 1 < < Mh^. Then, given any arbitrary F^ G J^h, 
either there exist an G J-^^ and a fCj^ G 7/i„ such that Fk = F^^ U Kj^ or 
there exist a i^^^ G 7/i^ and an G such that F^ = i^^-^ U F^^ . 

When considering functions in H^{Th}: we use the usual average and jump 
operators, respectively, defined for w G H'^^Th] along an interior face Fk by 



\ ((^li^i)li^fe + i.W\K2)\Fk) 



The above definitions are also valid for vector-valued functions w G [i/*{7/i}]^", 
in which case it follows that [w] ■ Vk = (w^li^JlFfc ■ J^Ki + {w\K2)\Fk " ^K2- 



3.2 Upwind Galerkin Approximation of the Vlasov Equation 



Here we describe the UG scheme for the Vlasov equation in full generality for 
a 2n-dimensional (n = 1, 2 or 3) phase space with inflow boundary conditions 
and piecewise polynomials of arbitrary degree approximating /. The simpler 
derivation of the method for a two-dimensional phase space with periodic 
boundary conditions in x and a piecewise constant approximation to / is 



given explicitly in Section |3.5[ since this method was used for the numerics of 
Section |4j Note, for this simpler version the derivation does not require use of 
the set of mesh faces J^h- 

For a given time T > and data trio {(y,fo,fj), the Vlasov equation along 
with the corresponding initial and boundary conditions are 



dtf + a-Vf = X (0, T] , (19a) 

f{t = 0) = f, n, (19b) 
/ = /, x(0,T], (19c) 
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where a = {v,—E) G M^", i? G is assumed given, V = (Va.,V,;), and 
Q = X Q^. It is assumed that = ^^=i[0, ^i] and = [— K, Vc]"", where 
Xi, . . . , Xn, Vc > are fixed. Following (|3]), we also define (keeping the same 
notation without loss of generality) the computational infiow boundary T^, 
associated with the computational domain Q, by 



= {(x,v) edn \ a-u <0}, 
where u is the outward unit normal to dQ. Then, = dQ \ F^. 



(20) 



For domains K C M^", let (-, ■)k denote the L^iK) -inner product. To distin- 
guish integration over domains K C M^""^, we use the notation (-, ■)^. A weak 
formulation for (19a)-(19c) is derived by multiplying equation (19a) by an ar- 
bitrary test function w G H^(Th) and integrating by parts over an arbitrary 
Kj E Th- This yields 



m, w)k, -(/,«■ Vw)k, + {fw, a ■ UK,)dK, = , V t G (0, T] 



(21) 



with being the outward unit normal to dKj and w denoting the interior 
trace of Kj, i.e., w~{x,v) = lim^^o- ((2;,^^) + si'Kj)- 



Upon summing equation (21) over all Kj and weakly enforcing the inflow 
boundary condition, we get 



k=l 



= - {frW,a-Uk)F,. (22) 

Approximating / by a function fh, which may be discontinuous across the 
interior faces, requires the introduction of a numerical flux A standard 
technique is to replace / by its "upwind value" /" on the interior faces [21j, 
where along each interior face is deflned by the given advecting flow 
field a{x,v,t) according to 

r{v,v,t;a) = \im f {{x,v,t) + sa{x,v,t))) = 

flKi{x,v,t) , if a{x,v,t) ■ Uk > 0, 
/|ii-2(x,w,t) , if a{x,v,t) ■ Uk < 0. 



Consistency follows from condition (23), since f^{a) = f on Fk whenever 
/ is continuous across each interior face F^. We also note that /" depends 
nonlinearly on a as, in general, if ai and 0:2 are two fiow fields having different 
values on F^ and if g E H'^{Th} is discontinuous across F^, then g^{ai +0^) ^ 
g^{ai) + g^i^a^)- Replacing / on the interior faces in (22) by its upwind value 
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leads to 

j=i k=i Ffeer^ 

which is the UG scheme for the Vlasov equation. Finally, defining the bilinear 
operator A by 

Nh Ph 

j=i k=i Ffesr^ 

(23) 

and the corresponding linear operator £, depending on the inflow data //, by 
C{w;ajj) = - {fiU!,a- iyk)Fk, (24) 

yields a variational formulation for the semi-discrete problem of flnding the 
fh e C^([0,T],D''^'''"''(7^)) approximation to /, satisfying, 



(25) 



ifhix,V,0),Wh)K, = {fo^Wh) 



(26) 



for all w G D'"""''"(7h), with /o and // approximations of the initial data 
f{x,v,0) and the inflow boundary data on Tj, respectively. 



We note that (25) produces a flrst-order ODE system to be described below. 
Indeed, for each Ki = Ki^ x Ki^ G Th, let ip\, . . . , ipl^^ be a basis for Q^'^^Ki^) x 
Q^^{Ki^), where rtb = [r^ + 1)" x [r^ + 1)" and then extend the domain 
of these functions to Q by deflning each to be identically zero in Q \ Ki. 
If j3{t) = {I3l{t), . . . , ^^^(t), . . . , ,9f'^(t), . . . , PnH^)) denotes the unique vector 
such that the UG approximation satisfles 



Nh rib 
j=l m=l 



(27) 



rib 



then fh^^ = J2 I^Li^L- Inserting (|27|) into (|25|) yields 

m=l 



Nh rib Nh rib 

EE ^^init)i^i,,whh + EE Piit)A{^i,,wh;^ 

j=l m=l j=l m=l 

= C{wh]aJ,), ywHED^-^'^^irH). (28) 
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Finally, since {'?/'p}pi'j^f=i is a basis for D^"='^"(7/i), then (28) is equivalent to 



rib 



j&N(i) m=l 

= C{tp;-a,f,), Wze{l,...,Nh},^pe{l,...,n,}, (29) 



771=1 



where N{i) contains the indices of all neighboring elements of Ki 



Equation ( 29 ) is seen to generate an equivalent matrix system, where rih rows 
of the matrix are generated at a time by sequentially taking i to equal 1, . . . , A^^ 
and for each i sequentially taking p to equal 1, . . . , n?, in Eq. (29). This proce- 
dure results in the matrix ODE system 



(30) 



where Ai is a constant matrix and A^ (a) is the corresponding sparse matrix, 
both of which are of dimension ribNh x n^Nh-, and L(a, /^) is a vector of length 

Since the support of the functions t/'^, . . . , ip\^ is Ki,Wi G {1, . . . , N^}, it follows 
that Ai is a block-diagonal matrix, where each block is an ra^ x matrix. This 
means that the inverse of Ai is easily computed. Thus, the UG approximation 
fh is equivalently defined to be the unique solution to 



-A^'A,{a)f3{t) + A^'L{a,f,] 



(31) 



where the initial condition /3(0) is uniquely determined by (26). To solve 
this system in time, a conservative explicit time stepping method such as 
the Runge-Kutta method can be used. 



3.3 Interior Penalty Approximations of the Poisson Equation 



In order to make this manuscript self contained, we also discuss the interior 
penalty formulations for the Poisson equation that weakly enforces approxi- 
mate continuity across interior mesh faces and Dirichlet boundary regions. As 
already noted, the mesh Th^ used to discretize the Poisson equation must be 
the same as that for the spatial domain used for the Vlasov equation, but the 
polynomial degree for the Poisson equation need not be equal to the degree 
in X for fh- 

Hence, for a given i) source G G L^(f2a;), ii) boundary data $^ e L'^{dflx,D) 
in the portion of the boundary referred as the Dirichlet boundary dflx,D, and 
iii) V$ ■ u = (homogeneous Neumann) or periodic boundary conditions on 
dflx \ d^lx,D, the boundary complement of the Dirichlet region, then the more 
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general form of Poisson equation with a positive definite permittivity term a 
discretized using the symmetric interior penalty (SIPG) method [TUfTT] . in- 
complete interior penalty (IIPG) method [BSll64f28j . or non-symmetric interior 
penalty (NIPG) [59] method is 



-V, -(aV,*) = G 



X 1 



(32a) 
(32b) 



where a{x) is any positive-definite continuous function in {C^{Qx))^^ 



For K C M", n=l,2, or 3, recall (■, ■)k is the L^(i^')-inner product with inte- 
gration over K C M"~^, while the notation (■, ■)^ is for boundary integrals. In 
addition, we use the identity 



dKi. 



+ Y: (V.<f (33) 



where uk is the outward unit normal to dK~ . 

Thus, the corresponding schemes SIPG, IIPG, and NIPG are all derived by 
multiplying (32a) by an arbitrary test function 9 G iJ'^jTh^}, s > 1/2, inte- 
grating by parts on each Kj^ G Th^, and summing each of the resulting local 
equations. Whence, one obtains the non-symmetric variational formulation 
given by 

Y: (aV.<l>, VJ)k,^ - Y ( «V.<I> [6] + [aV.<l>] 6, v^^ 

- (aV.$-z/,^,^^)^,^ = (G,^)n., (34) 



where the identity (33) was used to represent the inner boundary integrals. 

In particular, since the true solution for a bounded right-hand-side has at least 
the regularity $ G H^{Vt^) n H'^{ThJ, it follows that [$] = and [V^$] = 
along every interior face Fj.^ [32] • Thus, In order to get a good approximation 
for a regular solution satisfying these two jump conditions, one adds 'interior' 
penalty terms that vanishes on the true solution These terms are of the 
form 

Ph, 

where the values selected for Cg for the SIPG, IIPG and NIPG methods are 
-1, and 1, respectively. Similarly, such penalization is also required for the 
Dirichlet boundary region dQx,D as follows. 
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Indeed, we add to the bilinear form the interior penalty terms 

EjA^^.mMF.., (36) 

which are now completed by weakly enforcing both approximate continuity 
across the interior mesh faces and Dirichlet boundary regions dQx.D, by in- 
cluding the penalization 

E 7r^2(^-^o,0)F,^, (37) 

where a > is an arbitrary penalty parameter that is usually set equal to 
unity. Note that homogeneous or periodic boundary conditions vanish on 
boundary terms of the corresponding bilinear structure. Consequently, they 
do not require the boundary penalization term. 

Usually the penalization terms are denoted by the following non-symmetric 
bilinear form 

^^(^'^) = E 71-^ m),,^ + E ($,0)^,^ . (38) 



Adding both penalty terms and (35) to the left-hand-side of (34) results in 
Y: (aV.$, VJ)k,^ - E ( aV.<l> [6] + c^aVJ [$], u^^ )f,^ 

+ E 7)r^(^'^K = (G>^K+ E 7r^,^^u.e)Fk^. (39) 

Equation (p9|) completely defines each of the three interior penalty schemes. 



Setting AcX^i(^) equal to the first four terms on the left-hand-side of (39), 
where dependence on the parameter from (35 ) is noted as a subscript, we let 
the bilinear operator i3c^($, 6) := Ac^($, 6) + Jo-($, 6) and the linear operator 
H{9]G,^^) be equal to the two terms on the right-hand-side of (39). Then 
the function G D^-^{Th^) is the corresponding interior penalty Galerkin 
approximation to the Poisson solution $, if 



(40) 



Note, is positive definite (or coercive) even though each of Ac^{^,9) and 
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Jcr($,^) separately are only positive semi-definite [28|37f59l63|l64j . In partic- 
ular, Bc^ generates an equivalent norm for the Hilbert space H^iTh), 

if the measure of the Dirichlet boundary > 0. In the case of periodic 

or homogeneous Neumann boundary conditions we need the compatibility 
condition /^^ /ign fdxdv = \flx\, for each connected component fix of the 
spatial domain. In fact, the approximation to the Vlasov equation is done 
with a conservative UPG scheme and high order Runge-Kutta schemes, and 
in particular, for the case of periodic boundary conditions the compatibility 
condition at the numerical level is satisfied as well. 



Finally it is possible to see that the bilinearity and positive definiteness of 
Bcs for either a portion of Dirichlet or full Neumann or periodic boundary 



conditions, implies that (40) is equivalent to a uniquely solvable matrix system 
described next. Let = {r^ + 1)" and fi = 
be the unique vector such that 



> PI ; 



E E 



(42) 



Upon substituting this representation into (40 ), we conclude that (40 ) is equiv- 
alent to 

Bfi = H, (43) 

where B is an {ni, + l)Nh^ x [nf, + l)Nh^ invertible sparse matrix and H{G, ^d) 
is a vector of length (nf, + l)Nh^. However, B is not block diagonal, as was the 
case for Ai in the Vlasov ODE system. Thus, if the spatial domain is two- or 
three-dimensional, then using an iterative solver is in general the most efficient 



means of computing the solution fi in (43). However, for the one-dimensional 



spatial domain, fi can be computed by using an LU-matrix decomposition 
algorithm to factor B. For convenience, we write fi = B~^H, even though in 
practice fj, might be computed using an iterative method. 



3.4 Discontinuous Galerkin Approximation of the Vlasov- Poisson System 



The UPG method for the Vlasov-Poisson system results from combining the 
UG approximation of the Vlasov equation together with the interior penalty 
approximation of the Poisson equation. Thus, the approximation fh{t) to the 



solution f{x,v,t) of the Vlasov-Poisson system (2a)-(2b), at time t, results 
from an iteration as follows. 

Let fh(t) be given, where //i(0) is the approximation of the initial distribu- 
tion function. Then, an approximation ah{t) = {v, —Eh{t)) to a{x,v,t) = 
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{v, —E{x,t)) is determined by computing the corresponding approximation 
to —Eh{t), using one of the interior penalty approximation schemes to com- 
pute an approximation to the potential via the formula /i = 
B'^H{G,^d), where G is defined to be 1 - / fh{t)dv. 

Next, the approximate local field is computed by taking the local spatial po- 
tential gradients {Eh)\Kj^ = —^x{,^h)\Kj^ on each Kj^, which implies that 
Ehit) = Eh{fh(t)) is discontinuous across the interior faces of Th^. Conse- 
quently, it follows that the approximate flow field atit) = ah{fhit)), or equiv- 
alently ah(t) = ah{f3(t)), where f3{t) is a well defined given function. 

Summarizing, for any given fh{t), first compute the approximate ah{t), and 
then the approximate solution fh{t) to the Vlasov-Poisson system by solv- 



ing the ODE system (31) with a{t) replaced by ah{f3{t)). This leads to the 
following definition: 

Definition 1 The semi-discrete function fh{(3{t)) G C\[0,T], D''- is 
the UPG approximation to Vlasov-Poisson solution f , if fh{(3(t)) := fh as 



defined in (27) where f3{t) satisfies the nonlinear system of ODEs 



[tt) aMt)) = (v, V,<l>Mt))) with fi{t) = B-'H{l3{t)) , 

Pit) =-A^'A,{ah{/3{t)))/3{t) + A^'L{aMt)), f\) , (44) 



for all t e (0, T] . 



System (44) can be solved for (3{t) using any explicit time-integrator, following 
a classical Gummel map type iteration, where at any given approximation 
Z?""^, the /3(t"'~^) step (ii) is first computed since it does not involve time 
variation. Thus, one obtains a;/t(/3(t*~^)), which allows for the calculation of 
/?"■, the approximation to by step (iii). In our simulations we use a 

conservative high order Runge-Kutta time integrator. 

This very same iteration scheme was previously proposed in [H] for the calcu- 
lation of Boltzmann-Poisson solvers for semiconductors, and related work by 
the same authors cited below. There the calculation of the Poisson equation 
is done by a LDG scheme. 

We also point out that an error estimate for this nonlinear scheme can be 
found in [37|, which we state here in a concise form. Let (/, $) be a solution 



pair of the Vlasov Poisson system (2b), with boundary and initial conditions 



as described above, potential ■) G H^{^lx) for ^l^ C M", and distribution 
function / e G^{[0,T], H'^'iTh)) for n C R^"", for both s,s > n. Also, let 
denote the interior faces associated with element in Q, with F^ denoting 
the corresponding one associated with the elements in the 
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was defined at the end of Section 13. 1[ 



Furtlier, recall Q^lK) is the space of polynomials on a set K of degree less than 



or equal to r (cf. (16) and (17)), and and the degrees in x-space and v- 
space, respectively. Let the parameter A be the first eigenvalue to the Poisson 
equation in fl^, be the distributional solution to the perturbed Poisson 
equation for a source term ph, the charge density (f- average) associated with 
fh, and let p^ and be defined by 



Px = min{ra; + 1, s} and /x^ = min{r^ + 1, s} . 



(45) 



Then, we obtained the following error estimate for the 2n- dimensional semi- 
discrete formulation of the Vlasov Poisson system in terms the difference of 
suitable norms of potentials ^ — ^h, fields E — Eh = — V$ + V$/i, and particle 
distribution functions f — fh- 



1$ 



-2 



rT Ph 



+ [ lib/. ■ M'" [f - A]llo,ro + [ II l«/. • M"'[f - h]\\lr, 



where a is the penalization parameter of (36) and (37) and the || ■ H^/pc 
was defined in the previous subsection at (41). In addition, for a sufficiently 



smooth potential ■) G H^(^lx) for il^ C M" and distribution function 
/(t, ■) G H'^^{Q) for Q C M^", where the order of smoothness is given by the 
parameters s and s, this estimate is optimal. 



We close this subsection by noting that very recently our iteration scheme 
was reproduced in pl2] with a different Poisson solver. These authors perform 
error estimates for quadratic basis functions that preserve energy, but do not 
preserve the positivity of /. Numerical simulations have yet to be performed 
for their scheme and the amount of degradation cause by the lack of positivity 
remains to be ascertained. 
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3.5 Two-Dimensional Phase Space with Piecewise Constant Approximation 



We end this section with a description of the simphfied scheme for a two- 
dimensional phase space using piecewise constant approximations to /. As 
noted above, this is a positivity preserving (monotone) scheme that was used 
for the plasma simulations presented in Section |4j Such a piecewise constant 
basis function scheme can be easily extended to higher dimensions. We point 
out that, in work currently under preparation [TH|T^ . we extend the positivity 
condition to higher order basis functions by new limiter techniques inspired by 
[20ll74ll75j . which are maximum principle preserving and can be applied to both 
the Vlasov-Poisson and Vlasov-Maxwell systems. It remains a challenge to find 
a proper scheme that would preserved positivity and higher order moments, 
like momentum and energy. In the future we hope to compare our approach 
with extensive existing work on Vlasov-Ma xwell system |15|68|67j . 

Here, the simplified spatial and velocity domains are Qx = [0, L] and Q^j = 
[—Vc,Vc], with mesh points = a;o < xi < . . . < xn^^-i < xn^^ = L and 
-Vc = vo < vi < ... < vn,,^-i < vn^^ = Vc, where Nh^,Nh^ e N. Then 
for jx = l,...,Nh, and j„ = l,...,iV/,„, take = {Kjjf^^^ and = 



{Kj^y^Jl^i by defining each spatial element Kj^ = [xj^-i,Xj^] of size hj^ = 
Xj^ — Xj^-i, and each velocity element Kj^ = [vj^-i,Vj^] of size hj^ = Vj^ — 
Vj^-i, respectively. A mesh Th = {Kj}^J^i of the phase space domain Q is now 
generated according to Kj = Kj^ x Kj^ , where the index j is defined by the 
element ordering j = (j„ - 1)A^,,^ + j^, for = 1, . . . , Nh^ and = 1, . . . , Nh^, 
so that Th contains a total of = Nh^N^^ elements. The corresponding 
piecewise basis function is then given by setting 6'*"'(x) = 1, for x G -fCj^, and 
6''"=(x) = 0, otherwise, for = 1, . . . , Nh^. Similar basis is also constructed for 
i„ = 1, . . .,Nh,. Then, iakmg^\x,v) = 6'*-(x)x*"(^^), for i = 1, . . . , A^,,, 
generates the approximating space D^'^iTh) = span{-?/'^, . . . ,'?/'^''}. 



The corresponding upwind function /" defined in (23) on dKi \ F is now 
written in the simpler form 



,uf . N ) / {x,v,t) , li a{x,v,t)-VK, > 0, 
r{x,v,t;a) = { (47) 

/+(x,f,t) , if a{x,v,t) ■ UK, < 0, 



for f^{x,v,t) = lims^Q± f {{x,v) + suxiji) and the outward unit normal to 
Kj denoted by z/i^^(x,f) is simply defined by (0, —1) for v = (0, 1) for 

V = Vi^, (1,0) for X = Xi^ and (-1,0) for x = Xi^_i. 
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Therefore, the corresponding lowest order UG scheme is 
Ph 

k=l Fk&o 



In particular, for the piecewise constant UG approximation, fh{x,v,t) = 
J2 ip-'{x, v) to /(x, V, t), clearly one obtains {fh)\K = l^-'it), which implies 



and the corresponding semi-discrete UG approximation fh = J2f=i (^^ if) '^^ 
is the unique function in C^([0, T], D°'°(7h)) satisfying the initial condition 
Mx, V, 0) = Jj,^ /o, e {1, . . . , Nh} and Vt e (0, T], and 



JdKi/an JdK.nro 
+ I {fhm)))i(y^KjS = 0, for^ = l,...,iV,. 

(49) 

This last identity is a linear ODE system for any given electric field where 
the integration along the interior faces dKi \ F in (49) ( i.e., dKi fl dVt = 0), 
is simply 

/ fa-UK.dS = / E{x,t) {f{x,Vi^-i,t) - f{x,Vi^,t)) dx 
JdKi/an J^rx-i 



+ vifixi^,v,t) - f{xi^-i,v,t)) dv (50) 
and the integrations along dKi fl Tq and dKi n Tj satisfy 

/ fha-VK,dS = I {fh)i a - UK.dS = . (51) 
JdK.nro JdKinVj 

4 Numerical Results 



In this section numerical results are presented for six examples chosen to test 
the accuracy and convergence of the proposed DG method. The examples 
chosen are typical for testing Vlasov-Poisson algorithms (see e.g. [I^), but we 
have also included some atypical, more extensive comparison to theory. Four 
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of the examples test the hnear dynamics and its associated fine structure (fil- 
amentation) in phase space, while two examine the nonlinear evolution. The 
linear results are presented in Section 4.1 in Section 4.1.1 the ability of the 
DG method to solve the advection equation, the Vlasov equation with the 
electric field set to zero, is considered for both Maxwellian and Lorentzian 
equilibria, while in Section 4.1.2 the method is applied to the Landau problem 
and numerical results are extensively compared with the theoretical results of 
linear Landau damping, also for both Maxwellian and Lorentzian equilibria. 
The nonlinear results are presented in Section 4^ in Section 4.2.1 nonlin- 
ear Landau damping is considered while in 4.2.2 the nonlinear two-stream 
instability problem is computed. 



For all examples, piecewise constants are used to approximate the distribution 
/, piecewise quadratic polynomials are used to approximate the potential $, 
and time is discretized using a conservative fourth-order Runge-Kutta method. 
For all but the first two linear advection examples, the NIPG penalty method 
is used to approximate the Poisson system, and the linear system that results 
from using the NIPG method is solved using an LU-decomposition algorithm. 

Throughout this section, it is assumed that the distribution function / has 
the form 

f{x,v,t) = f,,{v) + 6f{x,v,t), (52) 

and the initial and boundary conditions used in all of the examples are of the 
form 



Sf{x,v,0) = Acos{kx) feq{v) , (53a) 
Sf{0,v,t) = 6f{L,v,t), (53b) 
<l>(0,t) = $(L,t), (53c) 

for {x,v,t) G [0,L] X [-Vc,Vc] x (0,T), where K > 0, L > 0, and T > are 
given. The constant Vc is the cutoff velocity and is chosen large enough so that 
the values of / are negligibly small when \v\ = Vc- It follows that each example 
is completely determined by specifying the governing equations along with the 
parameters feq{v), A, k, L, Vc, T. For both linear and nonlinear dynamics 
the initial condition is denoted by fo{x, v) = f{x, v, 0) = feq + Sf{x, v, 0). 



4-1 Linear Results 



Both the linear advection example of Section 4.1.1 with the initial condition 
/o = feq + Sf{x,v,0), where 6f{x,v,0) is given by (53a), and the Landau 
damping example of Section 4.1.2, governed by (13a)-(13c), require the speci- 
fication of feq. For both examples, the two choices for f^q introduced in Section 
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[2] are considered: the Maxwellian equilibrium /m of (14) and the Lorentz equi- 
hbrium of (15). 



Because it is most common to consider the Maxwell equilibrium, we explicate 
here several reasons for considering the Lorentz equilibrium, which to our 
knowledge has not been numerically tested previously in the literature. 



(1) Naturally occurring plasmas are sometimes not Maxwellian but posses 
kappa distributions [391H7] that have power law tails in v. The Lorentz 
equilibrium is in a sense an extreme case of these in that it has 
decay at infinity, with the existence of the particle density but not kinetic 
energy. In any event, distributions with power law tails are of physical 
interest and thus worth studying in their own right. (See also e.g. [65|66] .) 

(2) Because of the slow decay in v, the effect of truncating the velocity domain 
is amplified and a greater velocity domain is needed. This makes the 
Lorentz equilibrium a more stringent test for a numerical algorithm. 

(3) The linear dynamics of Vlasov theory is dominated by phase mixing. 



(4) 



the mechanism that underlies Landau damping (cf. Section 4.1.2). For 
the advection problem, the Lorentz equilibrium gives decay of the form 
exp{—kt), as opposed to the Maxwell equilibrium that gives decay of 



the form exp (— A;^t^) (cf. Section 4.1.1), and this suggests it might be 
a better test for getting Landau damping right. In fact, the reason for 
this exponential decay is that linear advection with the Lorentz equilib- 
rium shares the same analytic structure as that of the Landau damping 



problem (cf. Section 4.1.2), while linear advection with the Maxwell equi- 



librium does not. The essence of Landau damping can be traced to the 
Riemann-Lebesgue lemma ^2] , which states that charge density integrals 
of the form 



lim / dv q(v)e 



ivt 







provided g G L^, i.e. / dv\g{v)\ < oo. The rate of this temporal decay 
depends on the nature of the function g{v). 

The underlying reason for exponential decay in both the advection 
problem with the Lorentz equilibrium and the Landau damping problem 
with either equilibrium, is that the function g{v) for these problems is 
analytic in a strip in the complex f-space, and the damping rate is deter- 
mined by the pole closest to the real v axis. Thus, the basic mechanism 
of Landau damping is tested in the simpler advection problem when / is 
given by (52) with 6f given by (53a) and feq being the Lorentz equilib- 
rium. 

With the Lorentz equilibrium one can use residue calculus to explicitly 
obtain expressions for the damping rates (see Appendix |A]). Although for 
the advection problem this is also true for the Maxwellian, this is not the 



case for the Landau damping problem of Section 4.1.2 
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4.1.1 Advection Results 



The advection equation, 



0. 



(54) 



is a natural test for assessing Vlasov algorithms because the Poisson equation 
is removed from the calculation and the focus is placed on the resolution of 
phase space. With a Maxwellian equilibrium this example has been treated in 
many works, for example in pT|33|l54|[56] . In [56] four standard Vlasov solvers 
are compared. 

After computing /, the long time behavior of the solution can be checked 
by comparing the computational results with the known theoretical damping 
behavior due to phase mixing. To this end the net charge density ptot is given 
by 



ptot{x,t) 



-A 



— oo 
oo 



dv f{x, V, t) 
dv /o(x — vt, v) 

dv cos [k{x — vt)] feqiy) 



(55) 



where the second equality follows because the solution to (54) is given by 
f{x,v,t) = /q(x — vt,v), and the third upon substitution of (53a). According 
to the Riemann-Lebesgue lemma, limt_!.oo Ptot = under mild requirements on 
/eg. Below we give explicit forms for the decay for the Maxwell and Lorentz 
equilibria. 



Although it is common to consider the linear advection problem for test- 
ing numerical algorithms, it is not so well-known that there is an intimate 
relationship between the solution of the advection problem and the actual 
Landau damping problem. In fact, there is a one-to-one correspondence be- 
tween solutions of the two. In [5T|52] an invertible linear integral transform, a 
generalization of the Hilbert transform, called the G-transform, was explicitly 



constructed that maps (13a) into the advection equation (54). Thus, given an 



initial condition for the advection equation, there exists an initial condition 



for (13a) that transforms into the same solution. For the Lorentz equilibrium 



one can use residue calculus to obtain explicit expressions. In particular, the 



G-transform of the Lorentz equilibrium of (15) is 



G[fL 



VT 1 + f ^ 



1 + 



1 (1-3^;^) 
P (l + t;2)2 



(56) 



where the procedure is done mode by mode and k is the mode number. This 
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means that a solution to the linear advection problem with the initial condition 

fo = Acos{kx)h (57) 

is equivalent to the solution of the linear Vlasov-Poisson system with the initial 
condition 

fo = Acos{kx)G[h]. (58) 

The difference in long time decay between the advection problem and that 
of Landau damping can be traced to poles that occur in the integral trans- 
form. One could use the integral transform to further test the veracity of a 
numerical algorithm by comparing the solutions of the advection and Lan- 
dau problems, but we will not do so here. However, given the understanding 
provided by the integral transform, it is quite natural to examine the advec- 
tion problem with initial conditions that are meromorphic in velocity like the 
Lorentz equilibrium. 



Linear advection with a Maxwell equilibrium 



For our first example (54) is solved for feg = fu given by (14). We choose 



A = 0.1, k = 0.5, L = 4:71, Vc = 5 and T = 40. For this particular case, it 



is easily shown by elementary methods that the net charge density of (55) is 
given by 

pt^t{x,t) = -A cos (A;x) e-'^'*'/^ (59) 
This implies that max^. \ptot{x,t)\ = 0.1 e~*^/^, since k = 0.5. 

To test the accuracy and convergence of the UG method, max^^ \Ptot{x,t)\ is 
computed numerically and the results are plotted in Fig. [T] The numerical 
results were generated using the five uniform meshes {Nh^,Nh^) = (500,400), 
(1000,800), (2000,1600), (4000,400), (8000,400), where Nh^ and Nh, denote 
the number of partitions of the x-axis and the u-axis, respectively. 

The first three meshes are such that ~ h^, whereas the fourth and fifth 
meshes are such that ~ and hr^ ~ h^/lG, respectively. The motivation 
for using the last two meshes comes from the fact that the problem being ap- 
proximated involves only advection in the x-direction. Hence, it is reasonable 
to assume that mesh refinements in x will improve the numerical accuracy as 
much as performing refinements in both x and v, provided the refinement level 
in V is sufficiently small so that the error is almost entirely due to the refine- 
ment level in x. Figure [T] clearly shows that the UG method is both accurate 
and numerically convergent under mesh refinements. 

Our results compare with those of [TT|33f54ll56] for early times where solutions 
are accurate, but unlike the others we do not obtain the later time recurrence 
that arises from periodicity of f^^x—vt, v) in its first argument and the velocity 
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mesh size used in evaluating the integral of (55 ). This is because of the fineness 



of our mesh and because of the specific dissipative properties of the DG method 
give monotonic error. With a mesh of N^^ = 400 the recurrence time is Tr = 
27r/c/Aw = 7rNhJ{2Vc) ^ 126. 

Linear advection with a Lorentz equilibrium 



In this second example, we consider the advection problem with the Lorentz 



equilibrium (15). Numerical results are given for four different values of the 



wavenumber k. We will revisit these cases in Section 4.1.2, where we consider 



the actual Landau damping problem with the electric field not taken to be 
zero. Specifically, here Eq. (54) is solved for f^g = /l = 7r~^/(f^ + l), A = 0.01, 
and the large value = 30 for each of the wavenumbers k = 1/8, 1/6, 1/4, 
and 1/2. The corresponding values for L and T for k =1/8, 1/6, 1/4, and 1/2 
are L = 16n, 127r, 87r and An and T = 75, 75, 50, and 50, respectively. The 
uniform mesh {N^, N^) = (1000, 2000) was employed in each of the four cases. 
For this case, it is easily shown using residue calculus that the charge density 



of ( 55 ) is given by 



Ptot{x,t) 



-A cos (kx) 



-kt 



(60) 



This implies that max^, \ptot{x,t)\ = 0.01 e 



'kt 



The computed results for each of the four wavenumbers k are shown in Fig. [2] 
along with the exact result of (60). From the figure it is seen that for early 



times the computations match the theoretical result (60). At late times the 



computations diverge and, as anticipated, it is more difficult to resolve cases 
with larger k, i.e. with finer spatial structure. 




Fig. 1. (Linear advection with Maxwell equilibrium) Plot of log^^Q (maxx |/3tot(x, t)|) = 
-(logioe)t2/8 + 1 vs. t: analytic solution (solid), {Nh^,NhJ = (500,400) (dot), 
(1000,800) (dash-dot), (2000,1600) (dash- dot- dot), (4000, 400) (shoH dash), (8000, 
400) (long dash). 
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Fig. 2. (Linear advection with Lorentz equilibrium) Plot of logxo(niaxx |/3tot(x, t)|) = 
-(logioe)kt + 2 vs. t: k=l/8 (top left), k=l/G (top right), k=l/A (bottom left) and 
k=l/2 (bottom right); analytic solution (dash), numerical solutions (solid). 

4- 1.2 Linear Landau Damping Results 

The next two examples test the ability of the method to reproduce results 
consistent with the theoretical results of linear Landau damping. We empha- 
size that it is not enough to merely show exponential decay, but to believe an 
algorithm one must compare the decay rates and the parametric dependence 
of the theoretical rates. As above, both Maxwellian and Lorentzian equilibria 
are considered. For the Maxwellian equilibrium, which were previously treated 
in |llf33|l36f54f 731176] , results are computed using four successively refined uni- 
form meshes in order to demonstrate that the numerical decay rates converge 
to the theoretical decay rate under mesh refinement. For the Lorentz equilib- 
rium, we will see that the UG method is robust in the sense that it produces 
the correct decay rates for different wavenumbers k. 

As noted in Section |2| Landau showed the electric field decays exponentially 
in the time-asymptotic limit (for more rigor see for the linear case and 
the recent nonlinear results of [53J , of which Ref . 37 is an early version of the 
present work). More specifically, if we write the frequency as u{k) = Ufi{k) + 
ij^k), where Unik) and 7(/c) are real- valued, then in the time-asymptotic limit 
Ek{uj,t) decays at a rate 7(/c) and oscillates at a frequency uJuik). 
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Besides the Maxwellian equilibrium of ( 14 ) , the Lorentzian equihbrium fi of 



(15) gives rise to exponential damping of the electric field. In this case, the 



decay rate 7(fc) for the n-th electric field mode with k = 27rn/L is given by 

l{k) = -k, (61) 

where 7 < implies damping, and the corresponding frequency of the electric 
field oscillations satisfies 

ujRik) = 1 . (62) 
The derivation of (61)-(62) is described in Appendix [A} It is important to note 



that formulae (61) and (62) are explicit, which directly results from using the 
Lorentz equilibrium, whereas, as noted above, the formula for 7(/c) and u)R{k) 
when a Maxwell equilibrium is used are implicitly defined |42] . 



Linear Landau damping with a Maxwell equilibrium 



We solved the linear system (13a)-(13c) with f^g = Jm = (27r) -^'^e 



A = 0.01, k = 0.5, L = 47r, K = 4.5 and T = 80. For this problem, the 
theoretical decay rate and frequency of oscillations to the third decimal digit 
are respectively equal to 7 = —0.153 and lor = 1.415 [TT]. Approximate 
solutions are computed using the four uniform meshes {Nh^, Nh^) = (250, 400), 
(500,800), (1000,1600) and (2000,1600). 

Phase-space contour plots and cross-sectional plots in v of the approximate 
solution fh for {Nh,,NhJ = (500, 800) and t = 0, t = 25, t = 50, and t = 75 are 
displayed in Fig. |3} These sequential plots show the increase in filamentation 
of fh as time elapses. 

The convergence of numerical decay rates of the dominant electric field Fourier 
mode is demonstrated in Fig. |4j The decay rate resulting for each mesh was 
computed by calculating the slope of the straight line plotted in the figure. In 
each case, the line was defined by the point occurring at the peak of the third 
oscillation and the point occurring at the peak of the ninth oscillation. A time 
step of At = 0.001 was used in order to ensure that the points defining each 
of the lines were actual computed data points rather than points that were 
determined using some interpolation of the data. Under mesh refinement, the 
numerical decay rate is seen to converge, up to three decimal-digit accuracy, to 
the theoretical decay rate of -0.153. In all four cases, the numerical frequency 
is observed to correspond to the theoretical frequency up to three decimal- 
digit accuracy. We also note that, upon refining the mesh, the decay of the 
dominant mode is sustained for longer times before leveling off. Our results 
compare favorably with those of previous works |llf33ll36f54f 731176] , which were 
obtained by various other methods. Because of the fineness of our mesh we 
were able to proceed to longer times than all but [TH] which achieved machine 
zero for small perturbations. 
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As with the advection results above, in contrast to all other works, we do not 
see recurrence. It is important to note that recurrence is in fact an indication 
of numerical error. Recurrence is a general phenomenon in finite-dimensional 
dynamical systems with time advancement maps that are measure preserving, 
one-one, onto, bicontinuous, and have a bounded phase space. Poincare proved 
recurrence for finite-dimensional Hamiltonian systems, although it can hold 
for other systems as well. The Vlasov-Poisson system is an infinite-dimensional 
Hamiltonian system [12], and there is no general recurrence theorem for such 
systems. 

Numerical truncation procedures generally are not Hamiltonian, and indeed 
we have shown this to be the case for the DG method used here. However, 
it has been shown that the DG method, although not Hamiltonian, does give 
a finte-dimensional system that preserves phase space volume |18] and has 
a bounded energy. Thus our semi-discrete system has all the ingredients for 
making a general estimate for the recurrence time for the distribution function 
in terms of the number of degrees of freedom, like that done by Boltzmann for 
gas dynamics. This results is a very long time for meshes with any significant 
resolution. Note, this approach differs from a result that is often quoted in nu- 
merical works that follows for the method of [11] ; viz. for a given Fourier mode 
and a mesh of size Av the advection equation was shown to have a recurrence 
time for the electric field of Tr = 2nk/Av. Many Vlasov algorithms see recur- 
rence at a value that is close to this advection value Tr. It is important to note 
that this does not mean that the distribution function is recurring in phase 
space on this time scale. In [18] we present detailed recurrence calculations for 
our DG algorithm with piecewise constant and higher order polynomials. 

In any event, recurrence is not physical and the same can be said for the non- 
recurrent flattening of the electric field decay rate that is seen in our results 
with the DG algorithm at late times. The lack of recurrence in our results 
is due to both the very fine mesh size as well as the dissipative nature of 
DG algorithms, which is monotonic in nature. However, from Fig. |4] (and also 
Fig. |5] below) we would argue that DG can do a good job. 

Linear Landau damping with a Lorentz equilibrium 

The ability of the UPG method to achieve accurate damping results across 
different wavenumbers is now investigated. As noted above, the Lorentz equi- 
librium is used in this example because explicit formula for the electric field 
damping rates and the frequencies of the damped oscillations are easily ob- 
tained (see Appendix |A]) . Moreover, this example also tests the ability of the 
method to produce accurate results for an equilibrium that has a much heavier 
tail in v than does the Maxwell equilibrium. The heavy (algebraic) tail of the 
Lorentz equilibrium leads to a faster rate of filamentation because there are 
more resonant electrons than for the Maxwell equilibrium. 
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The linear system (13a)-(13c) is solved with feq = fi = 7r~^/(t'^ + l), A = 0.01, 
and Vc = 30. for each of the wavenumbers k = 1/8, 1/6, 1/4, and 1/2. Note, 
Vc needs to be large because of the algebraic tail of //,. The corresponding 
values for L and T for k =1/8, 1/6, 1/4, and 1/2, are L = IGtc, 127r, 87r and 
Arc and T = 75, 75, 50, and 50, respectively. Uniform meshes were employed 
in all four cases, where in each case {N^,N^) = (1000,2000). 

In Fig. [5} log plots of ptot for the four cases are given. The observed damping 
for k = 1/2 lasts for a shorter time duration than does the damping for the 
other three wavenumbers, even though the mesh for k = 1/2 has the finest 
resolution in x. This result is due to the fact that the filamentation in the 
velocity direction develops more rapidly than it does in the other three cases. 



From (61), obtained in Appendix |Aj, it follows that for a given wavenumber k 
the fundamental mode of the electric field damps exponentially in time at a 
rate equal to ■y{k) = —k and the frequency for all of the damped oscillations is 
u}{k) = 1. Therefore, the theoretical damping rates corresponding to k = 1/8, 
1 /6, 1/4, and 1 /2 are 7(fc) = -1/8, -1/6, -1/4, and -1/2. In all of the four cases 
shown in Fig. [5} the numerical damping rates and frequencies of oscillation 
are equal to the theoretical values up to the first two decimal digits. 
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Fig. 3. (Linear Landau damping with Maxwell equilibrium) Contour plots {lejf) 
and cross-sectional plots (right), x = 2-ir, for 6f at t = 0, t = 25, t = 50, t = 75 
(descending order). 31 




t 



Fig. 4. (Linear Landau damping with Maxwell equilibrium) Time decay plots 

of fundamental mode (with arbitrary ordinate origin) under mesh refinement: 
{Nh,,NhJ =(250, 200) (top left), (500, 400) (top right), (1000, 800) (bottom left) 
and (2000, 1600) (bottom right). The numerical decay rate converges to the theoret- 
ical value of -0.153 to within three decimal-digit accuracy; similarly, the numerical 
oscillation frequency agrees with the theoretical frequency of ur = 1.415 to within 
three decimal digits. 



32 



Numerical decay rate = -0.1 244 
Numerical frequency = 1 .001 




Numerical decay rate = -0.1 668 
Numerical frequency = 1 .002 




Numerical decay rate = -0.251 6 
Numerical frequency = 0.996 






Numerical decay rate = -0.5098 
Numerical frequency = 0.984 



20 30 
t 



Fig. 5. (Linear Landau damping with Lorentz equilibrium) Time decay plots of 
fundamental modes (with arbitrary ordinate origin): k=l/8 (top left), k=l/6 (top 
right), k=l/A (bottom left) and k=l/2 (bottom right). 
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4-2 Nonlinear Results 



Two nonlinear calculations are performed, an example of strong nonlinear 
Landau damping and a version of the two-stream instability that we examine 
in greater detail. In the former case, in addition to early time damping, we 
expect to see motions on the bounce time scale due to particle trapping, while 
in the latter we expect to see trapping and an asymptotic approach to a BGK 
state. 



4-.2.1 Nonlinear Landau Damping 

Following |lll)26ll31ll33l)41|l54|l73ll76j we consider an initial condition near a 



Maxwellian and of the form given by (52) and (53a)-(53c) with feq = f. 



Tm 



(27r)-i/2e-« /2, A = 0.5, k = 0.5, L = An, K = 5.0 and a run time of T = 120. 
Unlike the case considered in Section 14.1.21 we evolve under the full non- 
linear Vlasov equation. Solutions are computed using a uniform mesh with 
(iV,,,iV,J = (1000,750). 

Because A has been increased to 0.5, higher modes are excited at very early 
times and the damping rate significantly exceeds the linear rate of 7 = —0.153. 
This is because energy leaves the first mode as the higher modes get excited; 
i.e., in addition to the phase mixing process there is an energy transfer process 
caused by the nonlinear interaction of the modes. Figure [6] shows the ampli- 
tudes of the first four modes as a function of time, with a form similar to 
previous calculations. These plots are obtained from our mesh data by using 
the following 'log Fourier Mode' function: 

log FMk{t) = logio [\/\EsAt)f + \EcAt)f/L , ) (63) 

with 



and 



EsAt) ■= £dx E{x, t) sin (^^^ (64) 

Ec,k{t) ■■= £dx E{x, t) cos (^^^ , (65) 

where k is the mode number sought. From Fig. |6] it is seen that mode-one 
reaches its minimum value at around t 15 and then all modes grow until 
they reach their maxima at t ~ 40, consistent with previous calculations. 
Using the maximum amplitude of mode-one, the bounce time is calculated to 
be Tg ^ 20 and this is also in agreement with previous results. 

Examination of Fig. [7] reveals that we obtain a damping rate for the first mode 
of about 7 = —0.287, a value consistent with the —0.281 obtained by [TT] and 
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—0.243 by [73], given the different ways authors have used to malce this kind 
of fit. Also, given that we have a finer mesh, some deviation could be due to 
our more precise coupling to the higher modes. Examination of Fig. [8] shows 



there is significant entropy [(11) with C = —/In/)] dissipation at short times, 
as also seen by [33], and this introduces some error. Also note, we have run to 
T = 120, which is significantly longer than previous calculations and a small 
decay in all four of the amplitudes is seen. This could be due to transfer to 
higher modes, cascading, or due to dissipation in the algorithm. Over the full 



length of the run the total energy H of (10) is c onserved to within a few 
percent, and in the later part of the run entropy is well conserved. This is 
improved in the next section, where we treat the two-stream instability, by 
decreasing the mesh size. 

In Fig. |9] we plot the spatial average of the distribution function. Like other 
authors, we obtain early plateau formation in the vicinity of the phase velocity 
of the wave, seen in panel (b), which broadens as the higher order modes 
are excited. At around t = 40, approximately the time of the first bounce 
maximum, significant smoothing takes place and the system settles into a 
nearly constant average state with a persistent electron hole. The smoothing 
at t ~ 40 can also be seen in the results of [TT|3B] . Finally, we note the 
presence a small periodic dimpling behavior at the maximum that persists for 
late times. 



35 





Fig. 6. (Nonlinear Landau damping with Maxwell equilibrium) Amplitudes of the 
first four modes versus time. Mesh size {Nx,Ny) = (1000,750). 
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Fig. 7. (Nonlinear Landau damping with Maxwell equilibrium) Initial damping of 
dominant mode followed by growth due to particle trapping. 



(a) Energy. (b) Entropy. 

Fig. 8. (Nonlinear Landau damping with Maxwell equilibrium) The total energy H 



of (10), kinetic plus electrostatic, and entropy, (11) with C = —/In/, versus time. 
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4-2.2 Nonlinear Two-Stream Instability 



In this subsection, we numerically study the long-time nonlinear evolution of 
the two-stream instability, a standard example that has been used for demon- 
strating the efficacy of a variety of numerical algorithms 
Like these previous calculations, we consider the equilibrium state 



frsiv) 



2tt 



(66) 



and apply our algorithm to the nonlinear system (12a)-(12c) with data of 



the form given by (52) and (53a)-(53c). We choose parameters values to be 
A = 0.05, k = 0.5, L = Att, Vc = 5, and T = 100, a uniform mesh defined by 
{N^, N^) = (1800, 1400), and a time step At = 0.00125. 

Qualitative Behavior 



Figures 10 and 11 show 3D and 2D contour plots of the distribution function 
/ in phase space at different instances of time. The initial data consists of two 
symmetric, counter-streaming, warm beams with the small sinusoidal pertur- 
bation superimposed, as described above. We obtain the qualitative behavior 
expected: the linear two-stream instability grows exponentially until nonlin- 
earity becomes important and trapping occurs, with an eventual long time 
asymptotic Bernstein-Greene-Kruskal (BGK) state j8j with an apparent elec- 
tron hole-like structure (cf. e.g. ||60J). As the nonlinear system evolves, linear 



modes grow and saturate as shown in Fig. 12, the phase space hole forms as 



a portion of the electron distribution becomes trapped and exhibits filamen- 
tation. Over time, the sharp contour variation associated with filamentation 
is smoothed out, a consequence of nonlinearity and the numerical dissipation 
of the UPG method. (We note, here for aesthetic reasons some smoothing 
was also done by the plotting routine.) Our numerical results are indicative 
of a very stable computational scheme. This is partly due to the fact that our 
DG method is monotone and mass conservative. The particle number error 
accumulation in time is due to the nonlinear coupling. 

Early Growth 



As noted above. Fig. [12] shows the time evolution of the first four modes, which 
have wavenumber values k = 0.5, k = 1, k = 1.5, and k = 2, from t = to 
t = T = 100. (Recall, our domain is of size Air.) At early times, < t < 20, the 



behavior illustrated in Fig. 12 is consistent with the results of previous authors 

0.5 dominates the other modes that 
18 with growth 



where the initialized mode-one with k 
are not initially excited. All modes reach a maximum at t 
rates on the order of a couple of tenths of Up which is typical of linear theory 
for this problem. During early times, noise in the system and the growth of 
mode-one excites the other modes as was the case for the nonlinear Landau 
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damping of Section 4.2.1 A comparison was made in the calculations of |29j 



using results from [36j, but we give some further details here. 

The plasma dispersion function from Vlasov linear theory is given by 

1 fill 

where u will be off the real axis as is known for the two-stream instability. We 



evaluate this e with the two-stream equilibrium of (66) in order to obtain the 
linear growth rates. In Appendix |A] we show that e can be rewritten in the 
form 

e{z, fc) = 1 - ^ [l - + 2Z{z){z - z^)] , (68) 

where z = u/{k\/2) and Z is the usual plasma dispersion function as defined 
in [34j- Figure [13] is a plot of the growth rate obtained from e = adapted to 
our domain with L = An (cf. [36j Fig. 3) confirming that our early growth is 
about right. 

Conservation Properties 

Figure [14] shows plots of the invariants of Section [2] as functions of time, while 
Fig. [15] depicts the relative error of the total energy, {H — Ho)/Hq, and simi- 



larly the enstrophy. The top left panel of Fig. 14 shows that the total particle 
number is conserved quite well, with an error no larger than 0.01% over the 
full 100 units of time of our simulations. Actually, the DG discretization and 
Runge-Kutta (RK) method used for time advancement perfectly conserves this 
quantity for the transport equation. However the nonlinear iteration generates 
a monotone in time error for the total particle number of order of 10~^ in 100 



time units. The top right panel of Fig. 14 shows the evolution of the total 
momentum, P. This quantity is exactly conserved because the DG method 
cannot break symmetry, and the same is true for the time advancement al- 
gorithm. In fact the method perfectly conserves all odd velocity moments of 
/. The middle left panel of Fig. [14] depicts the evolution of the enstrophy 
Casimir invariant, / dxdvp, a q uantity that appears to be seldom-monitored 

it is seen 



in Vlasov codes ([33] being an exception). From the inset of Fig. 15 
to be conserved to within about 10%, which we suspect is comparatively good, 
and can be considerably improved by refinement and an increase in polynomial 
order. Conservation of this quantity arises because the solution of the Vlasov 
equation is a rearrangement, a property that we discuss below in more detail. 
It is violated because of small scale error and the diffusive nature of numerical 
algorithms. Note, like total particle number iV, the error in the enstrophy is 
monotonic in time. The middle right and bottom left panels of Fig. [14] show 
the evolution and error of the kinetic and electrostatic energies (§ and (ig. 



the first and second terms of (10), respectively. Individually these quantitie 



s are not conserved by the Vlasov equation, as is clearly evident from the 
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figures, but their sum, i7, shown in the bottom left panel, is conserved with a 
relative error less than 4% over the entire course of the run. The oscillations 
in the kinetic and electrostatic energies, indicative of the trapping process, 
cancel upon addition and give an error that increases monotonically in time. 
In Fig. 15 the relative error of the enstrophy and total energy are depicted. In 
other algorithms the error in H is oscillatory in nature, typically with small 
temporal mean. Even though this mean error could be small, it is important 
to realize that the oscillations amount, in a sense, to successively reinitializing 
the code, and the cumulative error of the solution for such a process is hard 
to assess. 



In addition to the conserved quantities of (|8]), (j9]), (10), and (11), the Vlasov 
equation possesses topological conservation. As mentioned above, it is well- 
known that the solution of the Vlasov equation is a rearrangement This 
means that the solution at any time t can be obtained formally from its initial 
value as follows: 



/(X, V, t) = fo{Xo{x, V, t),Vo{x, V, t)) = f^OZQ. 



(69) 



where Zq = {xQ{x,v,t),vo{x,v,t)) is obtained upon inverting the solution of 
the characteristic equations z = z{zQ,t). Because the characteristic equations 
have Hamiltonian form they preserve volume (here area), and 



d{xo,vo) 
d{x, v) 



(70) 



A consequence of this is the family of Casimir invariants of (11 ), whose invari- 



ance follows directly upon effecting the coordinate change zq z and making 



use of (69) and (70): 



dxQ / dvoC{fQ{xQ,Vo) = dx dvC{f{x,v,t) 



(71) 



where we have suppressed limits of integration. 



Under continuity assumptions, the level sets of /o are topologically equivalent 
to the composition /o o = /• This is what is meant by topological conser- 
vation. It follows that the number and nature of extrema, the values of / at 
these extrema, and the kinds of separatrices connecting saddle points must 
correspond to those of /q. In addition, although not a topological property, 
the area between any two contours of / must also be conserved, a consequence 
of the area preservation property of the characteristics. Since DG is a weak 
formulation with lack of continuity, the extent to which level sets are actually 
topologically conserved remains to be seen. Casimir invariants such as enstro- 
phy and entropy may be conserved well and be consistent with rearrangement 
inequalities such as Jensen's inequality without the continuity assumption. 
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Because the initial condition has the form fo{x, v) = /^^(f )[1 + 0.05 cos(x/2)], 
the initial phase space has nearly imperceptible initial 'tears' in the first 
panels of Figs. 10 and 11, that should be consistent with asymptotic state 
at t = T = 100. Extrema of fo{x,v) occur at points {x,v): (x,0), for all 
X G {0,47r}, (0, ±v^), and (27r, ±-\/2). Thus there is a trough along the line 
V = 0. The points (0, ±-\/2) = (47r, ±v^) atop the beams are maxima strad- 
dled by thin separatrices entering and emerging from the saddle points lo- 
cated at (27r, ±v^). The following extremal values of /o remain fixed in time 
under Vlasov dynamics, although their locations can change: /o(0, 0) = 0, 
/o(27r,0) = 0, /o(0,±v^) = 0.308, and fo{2^T,±^/2) = 0.279. A glance at 
Fig. [To] shows some degradation of the value of / atop the ridge-like struc- 
tures, while the vanishing value of / at the points (0, 0) and {2n, 0) is rigor- 
ously maintained. A measure of error in an algorithm is the extent to which 
it fills in the hole, i.e. returns values of /o(27r,0) 7^ 0. 



In our simulations the value of / remains zero at the minimum (27r,0) and it 
remains zero to high accuracy at (0,0). The values of / at the extrema and 
the associated separatrices atop the beams are never very discernible, their 
existence due to only a few percent in the variation of /, and are washed 
out by numerical dissipation as the contours of the distribution function wrap 
around and 'trap particles' over the course of the simulation. However, the 
trough at = is structurally unstable and separatrices emerge form here 
and connect the saddle point at w = and x = = 47r that straddles the 
minimum at (27r,0), the center of the electron hole. The eventual boundaries 
that delineate the trapped and untrapped particle populations, as the potential 
saturates into what appears to be a final BGK-like state, is what we discuss 
next. 



BGK Saturation 



Examination of Figs. 10, 11, and 12 shows the initial linear growth phase. 



followed by a particle trapping phase, and eventually a strong indication of a 
saturated state with clean contours of /, resulting at least in part from the 
small scale averaging inherent in any algorithm. It is generally believed that 
this evolution saturates, in some weak sense, to a BGK mode, although no 
proof exists. To test this belief we check to see if the contours of / are aligned 
with contours of the particle energy, S = — ^{x,t), which is well-known 
to be the case for an equilibrium state of the Vlasov equation (e.g. [8]). 

To this end we plot fioo{£{x,v)) := /(x,f,100), where in S{x,v) = — 
^{x, 100), the particle energy at t = 100. (Note, we suppress the time variable.) 
Figure [16] clearly indicates that a saturated BGK state has been achieved. 
Here, in the left panel, we have plotted / versus £{x,v), at t = 100, for all 
9 million pairs {x, v) of the uniformly distributed mesh over our phase space. 
Observe that to within the thickness of the line, /loo appears as a graph over 
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£{x,v) and that this is true even for large values of £^ where /loo > is 
maintained. Green and red dots correspond to positive and negative values 
of the velocity, respectively, and there does not seem to be any systematic 
directional bias. Because the computation was done for piecewise constant 
element functions, it is known that DG ensures f{x,v,t) remains positive for 
all times a nd this is then the case for /loo. The right panel of Fig. 16 is a plot 
of the electrostatic potential as a function of position for several instances in 
time as indicated in the figure. The curve labeled hj t = 100 is taken to be 
the saturated electrostatic potential that was used in the calculation of /loo- 



In Figs. 18 and 19 high resolution details of /loo are depicted. In Fig. 18(a) 
a cusp is seen at the trapping boundary that occurs at £^ = 0. A magnified 
view of this is shown in Fig. 18(b) In the course of the evolution the presence 
of the electrostatic force produces a trapped particle population with £ < 
0. It is conceivable that trapping, like scooping ice cream, is a non-analytic 
process and that this cusp represents something real about the mathematical 
nature of the dynamics. However, it also may be attributed to the choice of 
the numerical scheme, yet discontinuities have previously been proposed at 
trapping boundaries: in [60] a discontinuity is used in making electron hole 
models and in |30] two kinds of discontinuities are proposed for adiabatic and 
sudden trapping. None of these discontinuities match the cusp like feature 
that we observe, but the traveling wav e state treated in [3QJ is not the same 
as that reached by our simulation and it is possible that an analysis similar to 
that of [30] could produce a cusp. Further studies of this effect are currently 
being investigated. Figures [18(c) and 18(d) examine /loo near the minimum 
value of £. Note the clear steepening of /loo as it approaches its minimum 
value. It is possible that there is a universal nature to both this steepening 
and to the cusp at the trapping boundary, but we leave an investigation of 
this to future work. Figures 18(e) and 18(f) examine /loo near its maximum 



value. The spread here as well as the other panels of Fig. 18 give a sense of 
how close the system has relaxed to an equilibrium state. 



Figure 19 examines the details of /loo for higher values of £. In Fig. |19(b) 



splitting is seen to occur at around £ = 1.3. In Figs. 19(c) and 19(d)| we see 
that this small splitting persists to large values of It is interesting to note 
that because the red and green dots are mixed within each band, the splitting 
is not an effect of velocity direction. The splitting is within the resolution of 
the code and is believed to be a real effect, one that seems to indicate that 
complete saturation has not occurred. 



The above tests provide strong indication that the code has relaxed to near a 
BGK state. As further evidence we test to see if the charge density associated 
with /loo is consistent with that for a final BGK state. To this end we first 
observe that /loo can be fit reasonably well by a model distribution function 
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of the following form: 



(72) 



Here $a/ is the maximum value of and because £ = — $ we see — $m 
is the smallest possible value of £. This and positivity of / imply £* > ^m- 



Obviously (72) does not account for the fine features discussed above, but it 



will be sufficient for our purpose. 

The quantities a, $m, £*, and (3 are fitting parameters that can be matched to 
the code output. We obtain a rough idea of the degree of self-consistency, i.e. 



the degree to which Eq. ( 72 ) produces the correct electrostatic potential when 



substituted into Poisson's equation, by using the following values extracted 
from the code output: = 1-06 and £* = 1.59. For convenience we choose 
13 = 1. More significant figures in a, viz. a = 0.1148, are required to ensure 



that the total net charge is zero. In Fig. 17 we see these choices provide a 
pretty good fit to the data. However, insufficient charge near the maximum 
is spread to higher values of £. In a forthcoming publication we will examine 
this sort of modeling in much greater detail. 

We note, there is a relationship between £m, the value at which / obtains its 
maximum, and the other parameters. That is, f'{£M) = implies 



£ 



M 



2 - Pi^M + £*) + v/4 + /32($M - £*y 
2/3 



or 



£* 



<I>M - 13£m^m + 2£m - (3£li 



(73) 



(74) 



(3£m + 13^ M - 1 

where the latter is useful when £m is extracted from data. In fact, the value 



of £* = 1.59 used above follows from £m = 0.71 (cf. Fig. 18). 
Given /gt, the charge density is given by 

d£ht{£) 



/CO rc 
dvht = l- / 



(75) 



-*^2(£ + $) 

where in the second equality the formula £ = — $ is used to change the 



integration variable from v to £. With the choice (72) this quantity can be 
explicitly integrated to give 



Pfit = 1 - 



\/27r a 



I - I3\^M - $)($ - £*) + f ($M - 2$ + r 



(76) 



Defining a pseudopotential V by dV/d^ = pat, and integrating Poisson's equa- 
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tion yields 



where Vq is a constant and 



~2 



V($) = Vo, 



(77) 



V($) = $ 



(7J 



27ra 



^7/2 



15 



/3'(<I>A/ -$)($- £*) + -/3(<I>M - 2$ + r ^ 



is the expression for the pseudopotentiaL Here E"^ /2 acts as a pseudo kinetic 
energy. Thus we can interpret the BGK state by comparison to a particle 
dropped in the potential V at zero kinetic energy. Such a fictitious particle 
returns to a state of zero kinetic energy as x traverses the spatial domain, 
which must be the case if there is zero net charge. 



In Fig. 20, E"^ is plotted against $ for our simulation data. In light of (77) and 
(78), we expect E"^ to be a graph over $ at time t = 100, and from the figure 



this indeed appears to be the case. Also, since at t = 0, $" = Acos(x/2), it is 
easy to see that E^{x,t = 0)/2 = A^{x,t = 0) - $(a;,t = 0)78. This follows 
from the choice of ground for $ and the absence of net charge which gives 
via Gauss' theorem E{x = Att) = E{x = 0). With our boundary conditions 
E{x = 0) = 0. Thus we expect the parabolic curve in Fig. 20 labeled by 
t = with maximum $ occurring $jv/ = 8^ = 8 x 0.05 = 0.40. However, it is 
remarkable that at intermediate times E"^ is also a graph over $. This can be 
shown to be true if $ maintains symmetry about x = 2tt. This suggests that 
$ can serve as a good spatial coordinate, an idea we will discuss further in a 
forthcoming publication. 



The data of Fig. 20 can be compared to the model of Eq. (77) with (78). 
Choosing Vq so that ii^^($ = 0) = and the parameter values of Fig. 17 
yields the plot of Fig. 21 Here we see a reasonable fit to the data at t = 100 
Note, £^^($Af) ~ and the maximum value of E"^ is a close fit. We note, c 



certain level of accuracy is needed in the parameters because E is a sensitive 



function of $. We have not optimized the fit, but the result of Fig. |20]is roughly 
what one might expect for the expansion of /loo with only three terms of a 
Laguerre series, and thus gives a fair indication of self-consistency which was 
our goal. As noted above, we will revisit this again in great detail in a future 
publication. 
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Fig. 9. (Nonlinear Landau damping with Maxwell equilibrium) Spatial average of 
the distribution function at the times indicated. 
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Fig. 10. (Nonlinear two-stream instability) 3D plots of the solution / at t = 0, 
t = 20,t = A0,t = 60, t = 80 and t = 100. 
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Fig. 12. (Nonlinear two-stream instability) Time series plots of the first four modes: 
k = 1/2 (top left), k = 1 (top right), A; = 3/2 (bottom left), and k = 2 (bottom 
right). 




Fig. 13. Plot of growth rate 7 vs. k for the equilibrium distribution function of (66). 
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Fig. 14. (Nonlinear two-stream instability) Temporal evolution of the total particle 
number (top left), momentum P (top right), enstrophy (middle left), and kinetic 
energy (middle right), electrostatic energy (bottom left), and the total energy H 
(bottom right). See Section [2] for definitions of these quantities. 
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Fig. 15. (Nonlinear two-stream instability) Relative error of the enstrophy (left) and 
the total energy H (right). 
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Fig. 16. (Nonlinear two-stream instability) Plot of the distribution function at time 
t = 100, /loO) as a function of the particle energy £ = — ^{x,t = 100). Green 
dots correspond to positive velocities and red dots to negative velocities, v (left). 
Right panel depicts the potential $(x,t) at the times indicated. ^{x,t = 100) was 
used in £ for the plot of /loo (right). 
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Fig. 17. (Nonlinear two-stream instability) Model distribution function /gt of 
Eq. ([72]) with = 1-06, £* = 1.59, /3 = 1, and a = 0.1148 (solid blue) com- 
pared to code results at t = 100 (red dots). 
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Fig. 18. (Nonlinear two-stream instability) Fine scale details of the distribution 
function /ioo('?)> depicting a cusp formation near the trapping boundary at £ = 
(top left and right). Near the minimum value of £, /loo is observed to greatly steepen 
(middle left and right). The maximum of /loo is achieved at about £ = 0.71 (bottom 
left and right). In all plots red dots correspond to negative velocities and green dots 
to positive velocities, v. No asymmetry in the sign of v is evident. 
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Fig. 19. (Nonlinear two-stream instability) Plot of /loo(^) near E = 1.3 indicating 
mild splitting (toj) left), with stronger splitting seen at f ~ 1.5 (top right). The 
splitting continues to larger values of £ (bottom left and right). Red dots corre- 
spond to negative velocities and green dots to positive velocities v, with no evident 
asymmetry in the sign of v. 
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Fig. 20. (Nonlinear two-stream instability) Plot of E'^ versus ^> for the code output 
at the times indicated. Note, E'^ is a graph over $ even for the unsaturated states, 
i.e. times t < 100. 




Fig. 21. (Nonlinear two-stream instability) Plot of E vs. as given by (77) with 



(78) for the model /fit of (72) with the parameters of Fig. 12, and code output (red), 
both at time t = 100. 
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5 Conclusions 



In this paper we liave demonstrated tlie viability of the DG method with up- 
winding for solving the Vlasov-Poisson system. The method was described in 
detail for general polynomial bases of elements. Several examples were com- 
puted to demonstrate the utility of the method using a piecewise uniform 
basis. A convergence study was performed for the simple advection problem, 
indicating the degree to which filamentation can be resolved, and high res- 
olution computations of the linearized Vlasov-Poisson system were also per- 
formed. For linearization about a Maxwellian equilibrium, computation results 
were seen to compare very well with theoretical calculations of Landau damp- 
ing. Computations for Lorentzian equilibria were also presented and the wave 
number dependence of the Landau damping rate was verified, apparently for 
the first time in a simulation. The problems of nonlinear Landau damping 
and the symmetric two-stream instability were then considered, and results 
were compared to previous calculations. In both cases, constants of motion 
were monitored and the error was seen to be monotonic, due to nonlinear 
coupling. Also, reasons for the lack of numerical error in the form of recur- 
rence in our results were discussed. High resolution features of the distribution 
function were displayed for the long time BGK state that was reached in the 
two-stream calculations. Comparisons between theory and code results were 
made, particularly for the end BGK state, in an attempt to provide a high 
standard for truth in numerical code work. 

There are many future directions suggested by this work, some of which are 
ongoing. In a couple of upcoming publications [T8|T7] the authors will report 
on computations using higher order polynomial bases with an improved tem- 
poral integrator and a limiter that maintains positivity of the distribution 
function, as well as a study of numerical recurrence for these schemes. Indeed, 
[18] contains a thorough study of dependence of recurrence times and recur- 
rence amplitudes on the mesh size in x, v and the time step At, for Vlasov 
linear advection, as well as dependence on the order, type, and choice of ba- 
sis functions for the DG scheme. It is quite remarkable, that for the lowest 
order polynomial space of piecewice constant functions, one can prove that re- 
currence occurs, but the recurrence amplitude will decay, thereby suggesting 
value in choosing lower order functions. This fact significantly improves the 
performance criteria originally developed in 

As noted in the Introduction, the DG method can be run easily on nonuniform 
meshes, and in work not reported here we have seen this to be the case. This is 
a first step toward producing an adaptive code. Similarly, parallel implemen- 
tation awaits. There are many physical applications within reach, such as the 
treatment of a plasma diode which requires physical boundary conditions and 
the inclusion of various collision operators of relevance to plasma dynamics. 
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Since the Vlasov-Poisson system serves as a test bed for more sophisticated 
kinetic theories, it is our opinion that the DG method proposed here is an 
attractive alternative for many plasma physics computations. 
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A Appendix: Dispersion Relations 



Here we analyze the dispersion relations for the Lorentzian and two-stream 



equilibria given by (15) and (66), respectively. With the assumption that 
/(x, V, t) = feq{v) + Sf{x, V, t) and Sf{x, v, t) ~ exp{ikx — iut) and the scaling 
of variables described in Section |2} the plasma 'dispersion relation' is given by 



e{k, Lo) 



1 



{v — Lo/k) 



dv 



(A.l) 



where k is assumed to be real and positive and for bounded systems k = 
27Tn/L with n an integer. Landau damping arises upon analytically continuing 
this expression into the lower half w-plane and thus deforming the contour 
of integration |42j|. For integration along the real axis, stable and unstable 
eigenmodes (and embedded modes if they exist) satisfy 



e{k, u) = . 



(A.2) 



and it is this quantity we wish to investigate for both discrete modes and 
Landau 'modes', the latter obtained by analytic continuation into the lower 



half w-plane. In the latter case, the solution a;(/c) of (A.2) characterizes the 
asymptotic-time behavior of mode Ek{u,t) of the electric field. Specifically, if 
oj = uji+i'y, where ur and 7 are real-valued, then 7 < is the time-asymptotic 
rate of decay of the mode Ek^cu, t) and ur is the frequency of oscillation. 
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For the Lorentz equilibrium distribution function of (15), 



TT {V^ + 1) 



(A.3) 



Upon defining u = u/k and substituting (A.3) into (A.l), we obtain 

2 



e{k, ku) 



irk"^ J-oo (f^ + — u) 



dv 



(A.4) 



which expresses e as a function u and k. Evaluation of this integral is a straight- 
forward application of Cauchy's theorem. Assuming u is in the upper half 
plane, poles exist ai v = i and v = u, with corresponding residues of the 



integrand being Ri = i/[4:{u 



and Ru = u/{v? + 1)^. Summing over the 



residues gives the dispersion function 

e{k, ku) = 1 — 



1 



k'^iu + i) 



Upon setting e 



0, since uj 



ku, we obtain 



±1 — ik = oor + ^7 . 



(A.5) 



(A.6) 



which demonstrates, contrary to our assumption, that there are no discrete 
modes in the upper half w-plane. This is consistent with the well-known result 
that equilibrium distribution functions that are monotonic in f ^ possess no 



discrete growing or damped modes. However, continuing (A.5) into the lower 
half plane gives Landau damping with 7(fc) = —k and \ujR{k)\ = 1 for the 
time asymptotic behavior. 



For the two-stream equilibrium of (66), (A.l) gives rise to instability, i.e., for 



this case there is in fact a root with u in the upper half plane. For compu- 
tational reasons it is convenient to write e in terms of the plasma Z-function 
which is related to both the Hilbert transform and the error function (see e.g. 



pi]). Upon inserting (66) and performing some manipulation, (A.l) can be 
written as 

e = i-l[j,(z) + J2{z)] , (A.7) 



k^ 



where 



JAz) 



w dw 



w — z 



J2(Z) 



> dw 



w — z 



(A.8) 



and z = u/\/2. With the standard definition of the plasma Z-function, 



Z{z) 



TT J-oo 



dw 



w — z 



2ie-'' r e-''dt 

J — oo 



(A.9) 
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where in the first expression ^(-z) > and the value of Z for 3 (2;) < is 
obtained by analytic continuation, while the second expression is valid for all 
complex z. The second expression is desirable for computations. Also, some- 
times it is convenient to use the derivative formula 

Z' = -2(1 + , (A.IO) 

which is valid for all z. After some more-or-less standard manipulations and 



making use of (A.IO), we obtain 

2 



1 



A;2 



\-2z^ ^2zZ{z) 1 



(A.ll) 



We numerically evaluated this expression and searched for its roots to obtain 



Fig. 13, Because our system has the size L = Air, we write e{k,uj) with k 
replaced by k/2. 
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